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We study unconventional condensations of two- component bosons in a harmonic trap subject to 
the 3D a • p-type spin-orbit (SO) couphng. The topology of condensate wavefunctions manifests in 
the quaternionic representation. The spatial distributions of the quaternionic phase exhibit 3D 
skyrmion configurations, while those of the S'^ spin orientation possess non-zero Hopf invariants. 
As increasing SO coupling strength, spin textures evolve from concentric distributions to lattice 
structures at weak interactions. Strong interactions change condensates into spin-polarized plane- 
wave states, or, superpositions of two plane- waves exhibiting helical spin spirals. 
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Quantum wavefunctions are generally complex- valued. 
However, the ground state many-body wavefunctions 
of single component bosons are usually positive-definite 
known as '^no-node^^ theorem [T]. This statement is valid 
under very general conditions: the non-frustrated kinetic 
energy (e.g. the Laplacian type), arbitrary single particle 
potential, and coordinate-dependent interactions. This 
is a strong constraint, which implies that time-reversal 
(TR) symmetry cannot be spontaneously broken in con- 
ventional Bose-Einstein condensations (BEC), such as su- 
perfluid ^He and most ground state BECs of ultra-cold 
alkali bosons [3]. 

It is exciting to seek unconventional BECs with 
complex- valued condensate wavefunctions which sponta- 
neously break TR symmetry [4]. The "no- node" the- 
orem does not apply to spinful bosons with spin-orbit 
(SO) coupling because the kinetic energy is not just a 
Laplacian but also linearly depends on momentum. It is 
predicted that such condensates can spontaneously de- 
velop half-quantum vortex coexisting with 2D skyrmion- 
type spin textures [5|. Experimentally, spin textures in 
SO coupled exciton condensations have recently been ob- 
served [6] . The progress of synthetic artificial gauge fields 
in ultracold atomic gases greatly stimulates the investiga- 
tion of the above exotic physics 0, IH . Extensive studies 
have been performed for bosons with 2D Rashba SO cou- 
pling which exhibit various spin structures arising from 
competitions among SO coupling, interaction, and con- 
fining trap energy 

Quaternions are an extension of complex numbers as 
the first discovered non-commutative division algebra, 
which has^provided a new formulation of quantum me- 
chanics jl7l-[l9|. Can we have unconventional BECs with 
non-trivial quaternionic condensate wavefunctions? Sim- 
ilarly to complex numbers whose phases span a unit cir- 
cle, the quaternionic phases span a unit three dimensional 
sphere . Spin distributions associated with quater- 
nionic wavefunctions are obtained through the 1st Hopf 
map S'^. In 3D, both condensate wavefunctions 

and spin distributions can be topological non-trivial due 



to the homotopy groups 7r3(5'^) 7r3(5'^) = Z [20|, UM- 
The winding number of is the 3D skyrmion 

number and that of the S'^ is the Hopf invariant. 

Since one quaternion can be viewed as a pair of 
complex- number, two-component spinor wavefunctions 
map to single-component quaternion- valued ones. This 
provides a convenient view point to study its topological 
properties. In this article, we present 3D condensation 
wavefunctions with non-trivial topological configurations 
of quaternionic phases, which are stabilized by the 3D 
SO coupling of the a • p-type. The experimental realiza- 
tion of this SO coupling has been proposed by the au- 
thors through atom-light interactions in a combined tri- 
pod and tetrapod level system [22] and also by Anderson 
et al. 0] The condensation wavefunctions exhibit topo- 
logically non- trivial configurations as 3D skyrmions, and 
spin density distributions are also non-trivial with non- 
zero Hopf invariants. Spatial distributions of the quater- 
nionic phase textures and spin textures are concentric at 
weak SO couplings, and evolve to lattice structures at 
strong SO couplings. 

We consider a two-component boson system with the 
3D spin-orbit (SO) coupling of the <j • p-type confined in 
a harmonic trap. The free part of the Hamiltonian is 
defined as 
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where a =t,i refer to two internal states of bosons; A 
is the SO coupling strength with the unit of velocity; 
uo is the trap frequency. At single-particle level, Eq. [1] 
satisfies the Kramer-type time-reversal (TR) symmetry 
of T = {—i(j2)C with the property of = —1. However, 
parity is broken by SO coupling. In the absence of the 
trap, the good quantum numbers for single particle states 
are the helicity eigenvalues ±1 of a-p7|p|. This results in 
two branches of dispersions e±(k) = -^{k — ksoY where 
hkso = "rnX. The lowest single particle energy states lie 
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in the sphere with radius kgo denoted as the SO sphere. 
It corresponds to an SO length scale Iso = ^/^so in real 
space. The harmonic trap has a natural length scale It = 

y^^-, and thus the dimensionless parameter a = Irkso 
describes the relative SO coupling strength. We assume 
the contact s-wave scattering interaction defined as 
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and allow two different interaction parameters, includ- 
ing the intra and inter-component ones defined as g^^ = 
9ii = and g^i = cg^ where c is a constant. 

In the previous study of 2D Rashba SO coupling with 
harmonic potential [H, [l5| , the single particle eigenstates 
are intuitively expressed in momentum representation: 
the low energy state lies around a ring in momentum 
space, and the harmonic potential becomes the planar 
rotor operator on this ring subject to a 7r-fiux, which 
quantizes the angular momentum jz to half integers. 
Similar picture also applies in 3D [HI, [23|. The low 
energy states are around the SO sphere. In the pro- 
jected low energy Hilbert space, the eigenvectors read 
i^-\-{k) = (cos ^,sin ^e*"^^)-^. The harmonic potential is 
again a rotor Hamiltonian on the SO sphere subject to 
the Berry gauge connection as Vtp = ^m(zV/c — ^/c)^ with 
the moment of inertial / = Mkk^^ and Mk = /{muo'^). 
Ak = z(?/^+(^)|V/c|'0+(AJ)) is the vector potential of the 
U{1) magnetic monopole, which quantizes the angular 
momentum j to half-integers. While the radial energy 
is still quantized in terms of hu, the angular energy dis- 
persion with respect to j is strongly suppressed at large 
values of a as 
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where is the radial quantum number. As further 
shown in Ref. [22], in the case a ^ 1, all the states with 
the same rir but different j and jz are nearly degenerate, 
thus can be viewed as one 3D Landau level with spherical 
symmetry but broken parity. If filled with fermions, the 
system belongs to the Z2-class of 3D strong topological 
insulators. 

Now we load the system with bosons. The interaction 
energy scale is defined as Ei^t = gNo/l^ where Nq is 
the total particle number in the condensate. The corre- 
sponding dimensionless parameter is /3 = Eint/hu. The 
Gross-Pit aevskii (GP) equation at the Hartree-Fock level 
(HF) is obtained as 

r n^v^ ...^ ^ 1 
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where n^^^{f) = A/'o|V^^,|(r)p are the particle densities of 
two components, respectively. Eq. S] is solved numeri- 
cally by using imaginary time evolution. 



We first consider the limit of weak SO coupling, say, 
a ^ 1 and weak interactions. In this case, the gap be- 
tween single particle energy levels is large. The conden- 
sate wavefunction ijj remains the same symmetry struc- 
ture as the single particle one over a wide range of in- 
teraction parameter /3, i.e., ip remains the eigenstates of 
j = ^ as confirmed numerically, can be represented as 
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where /(r) and ^(r) are real radial functions, '^jij^ 
are defined as the SO coupled spherical harmonics with 
j = /±l/2 with the parity eigenvalues (— )^ respectively. 
For example, Y^^ = (1,0)^, and Y^^ = 

(cos 6>, sin l9e^^)^. The TR partner of Eq. [5] is V^^-,=-i = 
Ttjjj^j^^i. The two terms in Eq. [5] are of opposite 
parity eigenvalues, mixed by the parity breaking a • p 
term. The coefficient of i of the F-" term is because 
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the matrix element {Y^ \a ■ p\Y^ ) is purely imaginary. 
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For the non-interacting case, the radial wavefunctions 
up to a Gaussian factor can be approximated by spher- 
ical Bessel functions as /(r) ^ jo{ksor)e~^ /^^^ and 
^(r) ^ ji {ksor)e~^ /^^^ , which correspond to s and p par- 
tial waves, respectively. Both of them exhibit the pitch 
value around kso- At r = 0, /(r) is at maximum and g{r) 
is 0. As r increases, roughly speaking, the zero points of 
/(r) corresponds to the extrema of g{r) and vise versa. 
Repulsive interactions expand the spatial distributions of 
/(r) and ^(r), but the above picture still holds qualita- 
tively. In other words, there is a ^ -phase shift between 
the oscillations of /(r) and g{r). 




FIG. 1: The distribution of S{r) in a) the xz-plane and in 
the horizontal planes with h) z = and c) z/It = The 
color scale shows the magnitude of out-plane component 
Sy in a) and Sz in b) and c). The parameter values are 
a = 1.5, c = 1, and /3 = 30, and the length unit in these 
and all the figures below is It- 

The topological structure of condensate wavefunction 
Eq. [5] manifests clearly in the quaternion representation. 
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The complex two-component vector = {i/j^^i/j^)^ can 
be mapped to a quaternion variable through ^ = ^o+^i^+ 
whereto KeV^^,^i ImV^;,^2 -ReV^;,^3 = 
lim/j^; i^j^ k are the imaginary units satisfying P = p = 
— —1, and the anti-commutation relation = — = 
/c. The TR transformation on ^ is just — j^. Eq. [5] can 
be expressed in the quaternionic exponential form as 



a(r,(^) = \i(T)\e^^^^^^^^ = |^|(cos7 + cJsin7), (6) 



where |^(r)| = [f{r) +^^(r)]2, cos7(r) = /(r)/|^(r)|, 
and sin7(r) = ^(^)/|^(^)|. oj{Cl) = sinOcoscj) i + 
sin ^ sin ^ j + cos^ /c, which is the imaginary unit along 
the direction of Ct^ and satisfies = —1. Ac- 

cording to the oscillating properties of /(r) and ^(r), 
7(r) spirals as r increases. At the n-th zero point of 
g{r) denoted r^, 7(^n) = ^tt (n > and we define 
ro = 0), while at the n-th zero point of /(r) denoted 
<,7«) = (n-^)^ (n>l). 

The spatial distribution of the quaternionic phase 
gu;(Q)7(r) (^gfi]^g(j Eq. [6l which lies on S^^ exhibits a 
topologically nontrivial mapping from to S^, i.e., a 3D 
multiple skyrmion configuration. This type of topological 
defects are non-singular which is different from the usual 
vortex in single component BEC. For a closed 3-manifold, 
the Pontryagin index 7r3(6'^) = Z. However, the coordi- 
nate space is the open R^. At large distance r ^ It^ 
|^(r)| decays exponentially. The quaternionic phase and 
the mapping are not well-defined. Nevertheless, in each 
concentric spherical shell with < r < r^+i denoted as 
Shn, 7(r) winds from nir to (n + l)7r, and uj{Q) covers 
all the directions, thus Shn contributes 1 to the winding 
number of e'^^^^^^^^ on S^. If the system size is truncated 
at the order of /t, the skyrmion number can be approx- 
imated at the order of Irkso = In the previously 
studied 2D case, the spin density distribution exhibit 2D 
skyrmion configuration due to 712(8'^) [sl, IH, while 
condensation wavefunctions have no well-defined topol- 
ogy due to 7r2(5'^) = 0. 

Exotic spin textures in s pino r condensates have been 
extensively investigated [25-'275. In our case, the 3D spin 
density distributions S{f) exhibit a novel configuration 
with non-trivial Hopf invariants. S{r) can be obtained 
from ^(r) through the 1st Hopf map defined as S{r) = 
^i^t^a(3i^(3j or in the quaternionic representation. 
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where ^ = £,o—£,ii—£,2j—£,3k is the quaternionic conjugate 
of ^. The homotopy group from to 5^ is 7r3(5'^) = Z 
characterized by the integer valued Hopf invariant. The 
Hopf invariant of the 1st Hopf map is just 1 [21]. The real 
space concentric spheric shell Shn maps to the quater- 
nionic phase 5^, and the latter further maps to the S'^ 
Bloch sphere through the 1st Hopf map. The winding 



number of the first map is 1, and the Hopf invariant of 
the second map is also 1, thus the Hopf invariant of the 
Shn to 5*^ is 1. Rigorously speaking, the magnitude of 
S{f) decays exponentially at r /t, and thus the total 
Hopf invariant is not well-defined in the open R^ space. 
Again, if we truncate the system size at /t, the Hopf 
invariant is approximately at the order of a. 

Next we present numeric results for the spin textures 
associated with the condensation wavefunction Eq. [5] as 
plotted in Fig. [H Explicitly, S{f) is expressed as 
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In the X2:-plane, the in-plane components Sx and Sz form 
a vortex in the half plane of x > and 5*^ is prominent 
in the core. The contribution at large distance is ne- 
glected, where S{f) decays exponentially. Due to the 
axial symmetry of S{r) in Eq. [51 the 3D distribution is 
just a rotation of that in Fig. [T]a) around the z-axis. In 
the x?/-plane, spin distribution exhibits a 2D skyrmion 
pattern, whose in-plane components are along the tan- 
gential direction. As the horizontal cross-section shifted 
along the z-axis, S{f) remains 2D skyrmion- like, but its 
in-plane components are twisted around the z-axis. The 
spin configuration at z = — can be obtained by a com- 
bined operation of TR and rotation around the ?/-axis 
180°, thus its in-plane components are twisted in an op- 
posite way compared to those z = zq. Combining 
the configurations on the vertical and horizontal cross 
sections, we complete the 3D distribution of S{f) with 
non-zero Hopf invariant. 
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FIG. 2: The distribution of S{r) in horizontal cross- 
sections with a) z/It = —0.5, b) z//t = 0, respectively. 
The color scale shows the value of Sz, and parameter 
values are a = 4, /3 = 2, and c = 1. 

Next we consider the case of intermediate SO coupling 
strength, which suppresses the gap between the lowest 
and the second single particle levels. Interactions mix 
single-particle eigenstates with different values of j, such 
that rotational symmetry is broken and complex patterns 
appear. In this case, the topology of condensate wave- 
functions is still 3D skyrmion-like mapping from R^ to 
5^, and spin textures with the non-trivial Hopf invariant 
are obtained through the 1st Hopf map. Compared to 
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the weak SO coupling case, quaternionic phase skyrmions 
and spin textures are no longer concentric, but split to 
a multi-centered pattern. The numeric results of S{r) 
are plotted in Fig. [2] for horizontal cross-sections. In the 
x?/-plane, S exhibits the 2D skyrmion pattern: the in- 
plane components form two vortices and one anti- vortex, 
while Sz 's inside the vortex and anti- vortex cores are op- 
posite in direction, thus they contribute to the skyrmion 
number with the same sign. The spin configuration at 
z = zo > is shown in Fig. [2] (b), which is twisted 
around the z-axis clock-wise. After performing the com- 
bined TR and rotation around the y-axis 180 °, we arrive 
at the configuration at z = — zq. 



(a) 





FIG. 3: The distribution of S{r) in horizontal cross- 
sections with (a) z/It = 0, (b) z/It = 0.2, respectively. 
The color scale shows the value of Sz and parameter val- 
ues are a = 15, P = 0.8, and c = 1. The overall lattice 
exhibits the bcc structure. 

Next we consider the case of strong SO coupling, i.e., 
a 1. The single-particle eigenstates with = 
are nearly degenerate i.e., they are the lowest Landau 
level (LLL) states. We assume interactions weak enough 
to neglect the inter-Landau level mixing but still rela- 
tively strong to mix different states inside the LLL. Lat- 
tice structures are formed, and each site is a single 3D 
skyrmion for the condensate wavefunction ^(r), which 
corresponds to a spin texture with approximately a unit 
Hopf invariant. The numeric results for the spin texture 
configuration are depicted in Fig. [3] a) and b). In the 
xy-plane, the spin textures from a square lattice, and 
the lattice constant d is related to the SO length scale 
d 27rlso = 2t:It/(^' For the horizontal cross-section 
z ^ (i/2, spin textures also form a square lattice whose 
sites sit above the plaquette centers of the former lat- 
tice. Similar structure also occurs at z —d/2. Thus 
the topological defects form a body-centered cubic (bcc) 
lattice even though its size is finite confined by the trap. 

If in the strong SO coupling case, interactions are also 
strong to mix different Landau levels, we can neglect the 
effect of the harmonic trapping potential. The conden- 
sate wavefunctions select plane-wave states located on 
the SO sphere and are no longer topological. At c = 1, 
the interaction is spin-independent, and bosons select a 
superposition of a pair of states ±A? on the SO sphere, 
say, i/cgoi. The condensate wavefunction is written as 



with Na ^ Nb = No. The density of Eq. [9] in real 
space is uniform to minimize the interaction energy at 
the Hartree-Fock level. However, all the different parti- 
tions of Na^b yield the same Hartree-Fock energy. The 
quantum zero point energy from the Bogoliubov modes 
removes this accidental degeneracy through the "order 
from disorder" mechanism, which selects the equal par- 
tition Na = Nf). The calculation is in parallel to that of 
the 2D Rashba case performed in Ref. [HI, thus will not 
be presented here. In this case, the condensate is a spin 
helix propagates along z-axis and spin spirals in the xy- 
plane. At c 7^ 1, the spin-dependent part of the interac- 
tion can be written as Ha 



At c > 1, the interaction energy at Hartree-Fock level 
is minimized for the condensate wavefunction of a plane 
wave state e*^^°^| t), or, its TR partner. For c < 1, (Hsp) 
is minimized if (Sz) = in space. At the HF level, the 
condensate can either be a plane-wave state in the equa- 
tor of the SO sphere with spin polarizing in the x?/-plane, 
or, the spin spiral state of Eq. [9] with Na = N^. An "or- 
der from disorder" analysis on the Bogoliubov zero point 
energies indicates that the spin spiral state is selected. 
We also present the numerical results for Eq. |4] with a 
harmonic trap in Fig. 2] for the case of c < 1. The con- 
densate momenta of two spin components have opposite 
sign, thus the trap inhomogeneity already prefers the spin 
spiral state Eq. [9] at the HF level. 
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FIG. 4: The density profile (a) for ^-component, and 
that for ^-component is the same. Phase profiles for (b) 
t and (c) I- components, respectively. Parameter values 
are a = 10, ^ = 50, and c = 0.5. 

To summarize, we have investigated the two- 
component unconventional BECs driven by the 3D 
SO coupling. In the quaternionic representation, the 
quaternionic phase distributions exhibit non-trivial 3D 
skyrmion configurations from to S^. The spin orienta- 
tion distributions exhibit texture configurations charac- 
terized by non-zero Hopf invariants from to S'^ . These 
two topological structures are connected through the 1st 
Hopf map from to S'^. At large SO coupling strength, 
the crystalline order of spin textures, or, wavefunction 
skyrmions, are formed, which can be viewed as a gener- 
alization of the Abrikosov lattice in 3D. 

Note added Near the completion of this manuscript, 
we became aware of a related work by Kawakami et al. 
[28], in which the condensate wavefunction in the weak 
SO coupling case was studied. 
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